program schwarz_additif
  !
  ! Resolution du placien par l'algorithme de schwarz-additif
  !   ∆U = f
  ! 
  USE UTILS
  implicit none
  include 'mpif.h'

  ! Declaration des variables pour schwarz
  real*8, parameter  :: eps = 1.d-12
  real*8	:: ERR,ERR_glob
  
  
  ! Decalaration des variables MPI et charges	
  
  integer :: ierr, myrank, MPIsize
  integer :: linesByProc, remainder
  integer :: ibottom, itop
  integer status2(MPI_STATUS_SIZE)
  
  ! Declaration des variables numeriques
  integer	:: l,Max_l
  integer, dimension(:), allocatable :: nb_elem_proc
  real*8	:: t_debut,t_fin,t_max
  
  ! Declaration des variables du pb a resoudre
  integer	:: Nx,Ny,N, R
  real*8	:: dx,dy,Lx,Ly,posx,posy
  real*8	:: Aii,Cx,Cy,D
  real*8, dimension(:), allocatable :: U, Uold, RHS, Uschwarz, Urecv
  real*8, dimension(:), allocatable :: Utop, Ubottom
  integer	:: i,j,k
	

  call MPI_Init(ierr)
  
    
  call MPI_Comm_rank(MPI_COMM_WORLD, myrank, ierr)
  call MPI_Comm_size(MPI_COMM_WORLD, MPIsize, ierr)
  
  open(11, file ='param', status ='old')
  read(11, *) Ny
  read(11, *) Nx
  read(11, *) R
  read(11, *) Lx
  read(11, *) Ly
  read(11, *) D

  close(11)

  Ny = 10!*Np
  Nx = 10 
  !       R  = 1 !Recouvrement en Nbr de lignes (1 est le minimum)
  Lx = 1.0d0
  Ly = 1.0d0 
  D  = 1.0d0
  dx = Lx/(1+Nx)
  dy = Ly/(1+Ny)
  Aii = 2.d0*D/(dx*dx)+2.d0*D/(dy*dy) ! Terme diagonal de la matrice
  Cx  = -1.d0*D/(dx*dx)	! Terme extra-diagonal proche
  Cy  = -1.d0*D/(dy*dy)	! Terme extra-diagonal eloigne
  N  = Nx*Ny

  ibottom = 1
  itop = N
  
	! Initialisation du probleme
  U(1:N)    = 0.0d0
  Utop	    = 0.0d0
  Ubottom	    = 0.0d0

  linesByProc = Ny / MPIsize
  remainder = mod(Ny, MPIsize)

	if( myrank == 0 ) then	
		print*, "5\n"
	endif
	

	ALLOCATE(nb_elem_proc(1:MPIsize))


	if( myrank == 0 ) then	
		print*, "1\n"
	endif
	
	
	if(myrank == 0) then
		do i=0, MPIsize -1
			
			if (i < remainder) then
				ibottom = Nx * (i * (linesByProc + 1)) + 1
				itop = ibottom + Nx * (linesByProc + 1)
			else
				ibottom = Nx * (i * linesByProc + remainder) +1
		!~ 		ibottom = Nx * (i * linesByProc + myrank) +1
				itop = ibottom + Nx*linesByProc	
			endif
			
			nb_elem_proc(i+1) = itop - ibottom
			
		enddo
	endif


	if( myrank == 0 ) then
		print*, "2\n"
	endif
	
	
!~  Calcul des indices pour lme proc courant
  if (myrank < remainder) then
		ibottom = Nx * (myrank * (linesByProc + 1)) + 1
    itop = ibottom + Nx * (linesByProc + 1)
  else
		ibottom = Nx * (myrank * linesByProc + remainder) +1
!~ 		ibottom = Nx * (myrank * linesByProc + myrank) +1
		itop = ibottom + Nx*linesByProc	
  endif
  
  if( myrank == 0 ) then
		print*, "3\n"
  endif
  
  ! Allocation dynamique pour chaque proc en tenant compte du recouvrement
  ALLOCATE(U(1:N));ALLOCATE(RHS(1:N));ALLOCATE(Uold(1:N));
  ! Allocation dynamique pour la communication
  ALLOCATE(Utop(1:Nx));ALLOCATE(Ubottom(1:Nx))
  
  ALLOCATE(Urecv(1:nb_elem_proc(myrank + 1)))
  
  
  if( myrank == 0 ) then
		print*, "4\n"
  endif
  
	if( myrank == 0 ) then
		do j = 1, MPIsize - 1
			call MPI_Send(U(ibottom),nb_elem_proc(j+1), MPI_DOUBLE_PRECISION, j, 15, MPI_COMM_WORLD, ierr);
		enddo
	else
			call MPI_Recv(Urecv(1), itop - ibottom, MPI_DOUBLE_PRECISION, 0, 15, MPI_COMM_WORLD, status2, ierr);
	endif
    

  !	schw_iter = 0
  ERR_glob = 1.d3
  ! Remplissage du second membre incluant les conditions de Bords Dirichlet
  call second_membre(RHS,U,Utop,Ubottom,77,77,1,N,Nx,dx,dy,Cx,Cy)
  ! Resolution du systeme lineaire
  call jacobi(Aii,Cx,Cy,Nx,1,N,RHS,U,Uold,l)
  !call   grad(Aii,Cx,Cy,Nx,1,N,RHS,U,l)

  call wrisol( U,Nx,Ny,dx,dy,77,1,N )
  DEALLOCATE(U,RHS,Uold)
  DEALLOCATE(Utop,Ubottom)
  call MPI_Finalize(ierr)
  
end program schwarz_additif
